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Preface 


In  identifying  and  estimating  the  parameters  of  a  control  system  via 
quasilinearization,  the  output  "sentitivity  matrix"  is  perhaps  the  most 
costly  and  time-consuming  portion  of  the  required  calculation.  A  new  and 
highly  efficient  algorithm  for  calculating  the  sensitivity  matrix  of  a 
linear,  time  invariant  control  system  is  developed  in  this  paper.  It  is 
limited  to  the  single-input  single-output  case;  however,  it  may  easily  be 
modified  to  handle  several  inputs  and  outputs.  Also,  the  input  must  be 
piecewise  constant;  and  the  output  measurements  must  be  taken  at  constant 
time  intervals. 

Thanks  are  due  to  Dr.  J.  Gary  Reid,  who  has  researched  the  problem 
and  who  developed  the  basic  algorithm;  to  Dr.  David  A.  Lee,  who  helped 
me  with  the  mathematics;  and  to  Mrs.  Shirley  J.  Rapozo,  who  assisted  me 
with  the  typing.  And,  of  course,  special  thanks  are  due  to  my  husband, 
Leslie,  without  whose  computer  knowledge  and  moral  support  I  would  never 
have  finished  this  paper. 

Linda  K.  Palmer 
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Abstract 


In  this  paper,  a  new  algorithm  is  developed  to  calculate  the  output 
"sensitivity  matrix"  of  a  linear,  time-invariant,  single-input  single¬ 
output  control  system  with  piecewise  constant  input  and  output  measure¬ 
ments  taken  at  constant  time  intervals.  The  algorithm  incorporates  the 
singular  value  decomposition  to  investigate  parameter  identifiability  and 
estimation  accuracy  in  relation  to  the  system  as  a  whole  and  in  relation 
to  the  model  of  the  system.  As  a  result,  a  structural  condition  on  iden¬ 
tifiability  is  imposed;  and  the  system  designer  now  has  a  tool  to  evaluate 
how  well  the  model  describes  the  system. 

The  algorithm  is  verified  by  checking  its  results  with  those  using  a 
standard  software  package  for  numerical  integration.  It  is  then  used  to 
investigate  input  and  structural  design  issues.. 


I.  Introduction 


In  analyzing  a  system,  the  design  engineer  uses  a  mathematical  model. 
The  model,  by  its  very  definition,  represents  the  system.  It  must  account 
for  any  variations  in  the  system,  such  as  changes  in  the  initial  state  or 
in  the  input.  The  model  must  be  refined  for  the  uncertainties  in  the  sys¬ 
tem's  environment,  such  as  temperature  or  imperfect  measurements,  and  for 
the  effect  of  higher  order  terms  if  they  have  been  neglected.  The  parts 
of  the  system  affected  by  these  changes  and  uncertainties  are  called  the 
system  parameters . 

In  the  process  of  system  identification  and  parameter  estimation,  the 
model  parameters  are  found  and  given  values  (Ref  1).  This  paper  considers 
the  state  variable  model,  which  is  widely  used  since  it  requires  the  least 
amount  of  a  priori  information  to  predict  the  future  (Ref  15:2).  The 
state  variable  model  is  a  set  of  linear,  first-order  differential  equations. 
However,  the  relationship  between  the  system's  output  and  its  input  and 
model  parameters  is,  in  general,  highly  nonlinear.  (See  Eq  (4).)  There¬ 
fore,  system  identification  and  parameter  estimation  result  in  a  nonlinear 
optimization  problem,  where  convergence  is  desired  between  the  output 
predicted  by  the  model  and  the  actual  measured  output  of  the  system  (Ref 
15:3).  Basically,  the  parameter  identification  and  estimation  problem 
may  be  solved  by  linearizing  the  predicted  output  response  about  the 
current  best  estimate  of  the  parameters,  equating  the  result  to  the  actual 
measured  outputs  to  get  a  parameter  update,  and  repeating  the  process  un¬ 
til  the  parameter  update  goes  to  zero.  This  process  is  called  quasiline¬ 
arization.  For  more  details,  see  Eqs  (1)  through  (11). 
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The  most  costly  and  time-consuming  part  of  quasilinearization  is  the 
formation  of  the  output  "sensitivity  matrix."  (See  Eq  (8).)  This  paper 
is  concerned  with  the  formation  of  the  sensitivity  matrix  using  two  meth¬ 
ods.  The  first  method  involves  the  use  of  the  "sensitivity  system,"  a 
set  of  first-order  differential  equations,  and  is  treated  as  the  "stand¬ 
ard"  method.  In  the  second  method,  the  sensitivity  matrix  is  decomposed 
into  input-,  time-,  and  structure-dependent  parts,  and  may  be  referred  to 
as  a  "modal"  sensitivity  approach.  (See  Eqs  (33)  through  (59).)  The 
sensitivity  matrix  may  be  decomposed  because  the  system  output  may  be 
expressed  in  terms  of  the  system  input;  the  input,  output,  and  state  vec¬ 
tors  of  the  mathematical  model;  and  the  eigenvalues  and  eigenvectors  of 
the  system  plant  matrix.  The  reader  must  remember  that  the  formation  of 
the  sensitivity  matrix  constitutes  only  a  part  of  the  identification/ 
estimation  problem  and  must  be  calculated  at  every  iteration  of  the 
sol ution. 

In  Chapter  I.,  the  theory  of  the  quasilinearization  and  the  two 
algorithms  for  calculating  the  sensitivity  matrix  are  explained;  and  the 
algorithms  are  evaluated.  Solvability  of  the  problem,  or  i dent i fi abil i ty 
of  the  parameters,  and  the  accuracy  of  the  parameter  estimation  are 
discussed. 

Chapter  III  contains  the  computational  aspects  of  the  new  algorithm. 
This  includes  an  explanation  of  the  interactive  program  developed  to 
implement  the  algorithm  as  well  as  a  summary  of  the  computational  loading 
of  the  program.  Also  included  are  discussions  of  how  the  program  may  be 
used  by  the  designer  both  in  the  experimental  design  phase  and  in  the 
actual  estimation  task,  and  of  the  effect  of  variations  on  the  basic 
assumptions  of  the  algorithm. 
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The  new  algorithm  is  validated  by  comparing  its  result  with  that  of 
the  standard  sensitivity  system  method  in  Chapter  IV.  Then  the  new  algo¬ 
rithm  is  used  to  investigate  both  structural  and  input  design  issues  of 
the  model.  Several  examples  are  presented  and  analyzed. 

Conclusions  about  the  new  algorithm  and  recommendations  for  further 
research  are  made  in  Chapter  V. 
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II.  Theor 


Since  the  two  procedures  to  be  discussed  are  used  in  the  technique 
of  quasilinearization,  the  theory  of  this  technique  is  first  explained. 

The  singular  value  decomposition  and  the  insight  it  gives  into  identi¬ 
fiable  ity  and  accuracy  of  parameter  estimation  are  also  discussed. 

Then  the  "sensitivity  matrix"  S  of  the  state-variable  model  of  a  linear 
time-invariant  control  system  is  formed  first,  by  the  "sensitivity  system" 
method  and  second,  via  the  new  "modal"  method.  The  methods  are  compared, 
and  a  structural  condition  on  identi fiabil ity  is  presented. 

Quasi! inearization  (Ref  4:  6;  7 •  10;  11:  13:  14) 

The  parameters  of  the  system,  which  may  have  real  or  complex  values, 
are  assembled  into  a  vector  9_  of  dimension  NP.  The  single-input  single¬ 
output  system  considered  is  of  order  NA  and  has  state  dynamics 

x(t)  =  A(0)  x(t)  +  B(0)  U(t)  t  c  [0,  tf]  (1) 

x(0)  =  x  (0)  (?) 

—  — o  — 

with  output  measurements  at  each  time  t^ 
y(tk)  =  c(0)  x(tk) 

tk  £  to,  tf]  k  =  1  ,  2,  .  .  .  ,  K  (3) 

The  system  input  1)  is  assumed  to  be  known  exactly,  and  the 
nominal  values  of  A,  8_,  _C,  and  x^  are  assumed  to  be  good  approxima¬ 
tions  of  their  true  values.  The  matrix  A  is  assumed  to  be  non¬ 
defective  and  to  have  NA  distinct  eigenvalues.  All  system  matrices 
and  vectors  may  be  real  or  complex.  Parameters  e. ,  which  may  also 
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be  either  real  or  complex,  are  assumed  to  appear  linearly  in  the 
matrix  A  and  in  the  NA-dimensional  vectors  B,  C,  and  x^.  Each  parameter 
may  appear  more  than  once  in  the  plant  matrix  A  and/or  in  the  vectors 
B,  C,  and  x^. 

The  output  at  time  t^  may  be  calculated  from  the  equation 

At.  ft.  A(t .  -  t) 

y ( t k )  =  Ce  JQ  e  B  U(t)  dT 

NA  .  A.t. 

-  Z  (C  u,)  (v .  x  }  e  J 
S  =  1  J  * 


NA  .  A.t 

+  I  (C  u.)  (v.  B)  e  J 
j=l  J  J 


-A  .T 

e  J  U(t)  dT 


(4) 


where  u.  is  the  jth  eigenvector,  and  v.  is  the  jth  reciprocal  eigen- 
J  J 

vector.  (See,  for  example,  Reference  12,  page  3.)  These  vectors 
correspond  to  the  jth  eigenvalue  A.. 

J 

If  is  the  current  best  estimate  of  the  unknown  parameter  vector, 
then  the  predicted  output  response  may  be  written  as 


K®o> 


y(tr  y 

y(t2,  8^) 

y(tK’  V 

K  x  1 


(5) 
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+  H.O.T.  +  Noise  (6) 

Ignoring  higher-order  terms  and  noise,  Eq  (6)  may  be  written  in 
matrix-vector  form  as 

r  =  [  Y  -  YjeJ  ]  =  S  LB  (7) 

■>e  sensitivity  matrix  S  in  Eq  (7)  is  defined  by 

y(l  )(t1 )  y(2)(V  •  •  '  y(NP)(V 

y(l)^V  y(2)^V  '  ’  '  y(NP)^V 

S  = 

y(l)^K^  y(2)^K^  '  ’  ’  y(NP)^K^ 

where 

3y(tk) 

y(i)(V  =  30i 

for  k  =  l  ,  2.  .  .  .  ,K  and  i  =  1  ,  2,  .  .  NP  (Ref  15:  3). 

The  best  approximation  of  AS  is  obtained  by  finding  A?*,  the 
which  minimizes 

111  -  S  A6|| 2  00) 

If  S  has  rank  NP  (where  the  number  of  samples  K  is  assumed  to  be  greater 
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than  or  equal  to  the  number  of  parameters  UP),  then  a  unique  solution  £B* 
exists.  The  use  of  singular  value  decomposition  to  find  the  rank  of  S 
is  discussed  in  the  next  section  of  this  chapter.  The  calculation  of  A6* 
from  the  singular  value  decomposition  is  also  discussed  in  the  next 
section . 

Once  A6*  has  been  obtained,  the  parameter  vector  may  be  updated  by 

A  A 

B.  =6  +  AS * 

"new  ~°ol d  (11) 

The  process  is  repeated  until  convergence  is  obtained,  or  until  A0*  goes 
to  0. 

Singular  Value  Decomposition 

Earlier  in  this  chapter,  it  was  stated  that  S  must  have  rank  NP  for 
A9*  to  minimize  Eq  (10)(Ref  17:  634-635).  Readily  available  software 
packages  may  be  used  to  obtain  the  singular  value  decomposition  of  the 
matrix  S  (Ref  5:  LSVDF),  making  the  singular  value  decomposition  a 
reliable  and  easily  implemented  method  of  calculating  the  rank  of  S.  The 
singular  value  decomposition  also  gives  insight  into  the  accuracy  of  the 
parameter  estimation  and  into  the  'directions"  of  best  identification,  as 
will  be  explained  in  the  following  sections. 

The  Method.  If  the  matrix  S  has  rank  NP,  assuming  NP  £  K,  then  S 
will  have  NP  nonzero  singular  values  ,  .  .  .,  o^p  (Ref  17:  637). 

Then  S  will  have  the  singular  value  decomposition 
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where  E  is  the  diagonal  matrix  of  singular  values  such  that 

°1  -  °2  -  •  •  •  -  °NP  (13) 

(Ref  5:  LSVDF).  The  orthogonal  matrices  UT-j  and  V  are  composed  of  the 
left  the  right  singular  vectors  of  S,  respectively.  The  unique  vector 
A8*  is  then  given  by 

=  V  r1  UT1  r  (14) 

See  Stewart  (Ref  17;  636)  for  mathematical  details. 

Since  E  is  diagonal,  E~^  is  formed  by  taking  the  multiplicative 
inverses  of  each  of  the  elements  along  the  diagonal  of  E.  If  E  is  not 
invertible,  then  some  of  the  singular  values  of  S  are  identically  zero. 

Then  V~ ,  the  columns  of  V  corresponding  to  the  zero  singular  values, 
span  the  null  space  of  S,  and  the  unique  vector  A6*  minimizing  Eq  (10) 
does  not  exist.  This  means  that  not  all  of  the  parameters  of  Eqs  (1) 
through  (3)  can  be  identified,  and  the  system  is  termed  "noni dent i fiabl e" 
(Ref  15:  10).  When  this  occurs,  the  system  must  be  modified  either  by 
eliminating  or  combining  several  of  the  unknown  parameters  or  by  remodel¬ 
ling  the  system  itself  (Ref  15:  5-7)  so  that  S  has  rank  NP. 

Accuracy  of  the  Estimation.  Once  the  parameters  have  been  estimated, 
a  measure  of  the  quality  of  the  estimation  is  desired.  To  find  a  suitable 
measure,  Stewart  derives  the  following  spectral  norm  relations  (Ref  17:  636) 


=  sup 


°NP  "  Gmin 
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where  ||  '  ||  denotes  the  spectral  norm  and  S  is  the  pseudo-inverse  of 
S,  defined  by 

ST  S  S  =  S  (17) 

and 

S  S7  S  =  S  (18) 

(Ref  17:  634).  Then,  using  Eq  (12),  S  may  be  found  to  be  equal  to 

s"  =  V  r1  UT-j  (19) 

According  to  Stewart,  the  solution  to  the  perturbed  problem 

£  +  6  =  S  (A6  +  e)  (20) 

satisfies  the  error  bounds 


(21) 


where  T-|  is  part  of  the  partitioned  vector 

r  =  [T1  r2]  (22) 

and  where  k,  the  ratio  of  the  largest  to  the  smallest  singular  value  of 


S,  is  called  the  condition  number  of  S 


°2  (23) 

The  condition  number  is  bounded  by  unity  from  below  and  infinity  from 
above.  The  nonnegative  constant  o  is  defined  by 

II r-j!|  =  n  01  II A6  II  (24) 

and  is  problem-dependent.  It  may  be  shown  that 


— —  <_  r)  £  1 

K 


(25) 


When  S  is  ill-conditioned,  i.  e..  when  S  is  large,  then  n  may  be  close 
to  k”1  (and  both  are  near  zero).  Then  T  is  said  to  "reflect  the  ill- 
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condition  of  S."  If,  however,  n  is  near  unity,  then  r  does  not  reflect 
the  il 1 -condition  of  S.  Then  the  estimation  accuracy  is  proportional 
to  •'  or  to  some  power  of  k,  making  the  condition  number  a  convenient 
and  reliable  measure  of  the  condition  of  S  (Ref  17:  653-655).  No  matter 
what  n  is,  however,  k  will  always  give  a  worst-case  bound  on  the  error: 

The  closer  k  is  to  unity,  the  better  is  the  parameter  estimation.  When 
<  is  not  close  to  unity,  the  accuracy  of  the  estimation  may  be  question¬ 
able;  and  perhaps  the  model  should  be  changed. 

"Directions"  of  Best  Identification.  In  the  singular  value  decom- 

oosition  defined  by  Eq  (12),  the  NP  columns  of  V,  u . ,  j  =  l  ,  2,.  .  .  ,NP, 

NP 

rorr  an  orthogonal  basis  for  the  parameter  space  R  such  that 

!! S  Vj  II  2  =  u'  (s's)  v.  =  a2  (26) 

Since  the  singular  values  are  ordered  from  largest  to  smallest,  gives 
the  "direction"  of  best  identification;  while  gives  the  "direction" 
cf  worst  identification.  This  means  that  the  first  parameters,  corre- 
s  p : r  d i n  g  to  the  larger  singular  values,  will  be  easier  to  identify  and 
esvnate  accurately  than  the  last  parameters,  which  correspond  to  the 
r~a‘ler  singular  values.  If  the  estimate  is  too  inaccurate,  the  designer 
r i c h t  try  to  reduce  the  model,  thereby  reducing  the  number  of  unknown 
caraneters  ( Ref  8) . 

The  next  two  sections  of  this  chapter  are  concerned  with  calculating 
the  sensitivity  matrix. 

:er;itivity  System  Calculation  of  S 

Toe  ’standard"  sensitivity  system  method  of  calculating  the  sensitivity 
"at’-iy  is  described  in  this  section.  For  ease  of  notation,  functions  of 
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£,  e.  g.,  A=  A(0_),  are  implicit.  In  addition,  x_  and  U  are  functions  of  t. 
Taking  the  partial  derivatives  of  Eqs  (1)  and  (3)  with  respect  to  the  ith 
parameter  yields 


-(i) 

y(i)(V 


A  -0)  +  A(i)  -  +  -0)  U 

I  x(l-)(tk)  +  C(i)  x(tk) 


for 

i=  1, 

2,. 

X 

-0) 

-(NP) 

where  the 

new 

(27) 

(28) 

K.  The  sensitivity  system  becomes 


0) 


o 

A 


(NP) 


0 


.  0 

0  A 


X 

B 

-(1) 

-(1) 

• 

+ 

• 

1  J* 

-(NP) 

(29) 


where  the  new  "plant  matrix"  is  formed  by  placing  A  along  the  diagonal, 
the  partial  derivatives  of  A  along  the  first  column,  and  zero  matrices 


el  se 

where . 

Similarly,  the 

new  "output" 

equatio 

~y 

■C  0  • 

•  •  o  - 

-x 

y(D 

■ 

-(1)  -  ‘ 

•  •  0 

-(1 ) 

1  • 

j 

■ 

•  0 

i_y(:iP)_ 

o 

CL. 

c3T 

•  0  c 

J* 

“O 

1 _ _ 

(30) 


for  each  time  t^=  1,2,.  .  .,  K.  The  differential  equations  (See  Eq  (29).) 
nay  be  solved  using  ODE  (Ref  10)  or  a  similar  software  package. 

The  difficulty  with  this  method  is  that  it  results  in  at  most 
Vxr  *  (NP  +  1)  and  at  least  NA  *  (NPA  +  NPB  +  NPX  +  1  -  number  of  parameters 
cc'vnon  to  A,  B_  and  x^)  coupled,  linear  differential  equations  which  may  be 
"stiff"  and,  therefore,  hard  to  integrate  accurately.  This  occurs  when 
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the  eigenvalues  of  the  plant  matrix  A  are  far  afart.  An  example  of  a 
very  stiff  system  is  the  two-dimensional  case  '"A-  ?)  where  fe  eigen¬ 
values  A-j  and  equal  -1  and  -100,  respective ;  y. 

Another  shortcoming  of  this  method  is  its  "ur.contrc!  laf  i  1  i  ty . " 
occurring  when  the  number  of  equations  in  the  syste'.  is  too  large.  1 
Unnecessary  work  is  created  in  a  direct  integration  of  the  excess 
equations  of  the  system.  To  lessen  the  workload,  model  order  reduction 
techniques  may  be  used  to  reduce  the  number  of  equations  from  ’.A  *  (NP  1) 

to  about  2  NA  for  a  single-input  single-output  system  (or  to  about  2  ‘.'A 
*  MB,  where  NB  is  the  number  of  inputs  of  the  control  system)  (Rc-f  C). 
Unfortunately,  the  linear  transformations  required  by  the  model  order 
reduction  techniques  may  significantly  increase  the  error. 

In  addition,  this  method  gives  no  indication  about  how  well  the 
model  described  by  Eqs  (1)  through  (3)  actually  represents  the  control 
system.  The  desirability  of  such  "structural"  insight  leads  naturally 
to  the  decomposition  of  S  into  time-,  input-  and  structure-dependent 
parts . 


Modal  Algorithm  for  Calculation  of  S 

In  this  section,  a  new  algorithm  for  calculating  the  sensitivity 
matrix  is  developed.  Because  the  output  sensitivities  are  expressed  in 
terms  of  the  eigenvalues  and  eigenvectors  of  A,  this  method  may  be  called 
the  "modal"  approach: 

Frcm  linear  systems  theory,  when  the  NA-dimensional  matrix  A  has 
NA  distinct  eigenvalues 


ft  +■ 

At 

e 


k 


NA 


(31) 
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t 


Noting  Eq  (31)  and  taking  the  partial  derivatives  of  all  the  terms  in 
Eq  (4)  with  respect  to  0.  gives 


NA 


•Yk 


X.t. 
J  k 


hn'V  *  ^  [  %f>  h  e  h  +  £ih(i)  e 

M'k  -  Vk  - 

+  C  u .  e  J  v.  x  +  C  u .  e  v  .  x 

--J  -J(1)-o  --J  -J^(i) 

X.t  .  At  . 

+  C/ . s  u  .  e  k  v .  B  +  Cu.  eJ  v,B 
-( i )  -J  -J  -  -  “J  ( i )  “J  ~ 

X.t  -  X.t,  „ 

+  -  -J e  J 


x 

-o 


r  ' k  -X  .T 
J  e  3  U ( t )  di  ]  (32) 


Calculating  y(YY  ^or  i  =  1  >  2 , .  ■  NP  and  k  =  1,  2,.  .  K  gives 
the  elements  of  the  K  x  NP  sensitivity  matrix  defined  by  Eq  (8).  By 
rearranging  the  terms  on  the  right-hand  side  of  Eq  (32)  for  all  i  and  k, 
and  by  expressing  the  result  in  mat ri x- vector  form,  the  sensitivity  matrix 
may  be  written 

S  =  [E  F]  G  (33) 


where  E  depends  on  time  alone,  F  is  dependent  on  both  time  and  input,  and 
G  depends  solely  on  the  structure  of  the  system  matrix  and  vectors  (Ref  15: 
4).  The  matrices  E  and  F  have  row  vectors 


Yk  Yk 

e  ,  e 


'NAtk  .  oXltk 

e  *  t.  e 

k 


.  Yk  .  :NAtk 

tk  e  tk  e 
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1  x  2.NA  (34) 
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with 


FI  0.-0=  C4*Ud(1)  '  (44) 

F3( 1 ,1)  =  C5*Ud(l)  (45) 

A.A 

where  A  equals  the  sample  spacing,  C4  =  (e  -  1)/A„  .  and  C5  =  (1  +  AA« 

A.A  AfA  2  £  l 

*e  -  e  )/A^,  .  See  Appendix  A  for  the  derivation  of  Eqs  (38)  through 

(45).  Note  that  since  appears  in  the  denominators  of  constants  such 

as  C4  and  C5  (See  Subroutine  EFMAT  in  Appendix  B.),  all  the  eigenvalues  of 

A  must  be  nonzero. 

The  input-  and  time-dependent  matrix  [E  F]  has  dimension  K  x  4NA, 
where  K  is  the  number  of  samples  taken.  This  matrix  is  multiplied  in  Eq 
(33)  to  the  4NA  x  NP  structure-dependent  matrix 


G  = 


Gzi 


zs 


4NA  x  NP 


(46) 


where  zi  stands  for  zero  input  and  zs  stands  for  zero  state.  The  ith 


column  of  is  given  by 


V.<*i  *o>] 

seT^c  (^a 

1  x,)] 

if4A) 


f 
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[-C-(i)  y_i  +  C  u1(i)]  (y-!  +  (Cu,)  [V1(1J  ^  +  v1  ^(i)] 


(i) 


CCM\  “MB  +  3  (Vjv.fi  *J  +  (C  Uk,J  [v 


(i)  T  —  ^A/  .  v J  V^A  ±o>^±2HA>  l^NA/,»V  ^1A  V. 


( i ) 


(t: 


( i ) 


X,  (C  u, )  (v.  x  ) 
( i )  ~_1  _1 


>-  (C  u,M)  (v^  %) 


'2(0 


(47) 


The  ith  column  of  G  is  the  same  as  that  of  G  .  except  that  x  is 

ZS  Z1  K  -0 

replaced  by  E[  (Ref  15:  4,  5).  At  this  point  B_.  C,  x  .  u.  and  y.  are  known. 

P  J  J 

The  sensitivities  of  the  quantities  necessary  to  compute  G^  and  G2$  are 
discussed  in  the  next  two  sections. 

Eigenvalue  and  Eigenvector  Sensitivities.  Crossley  and  Porter  have 
derived  closed  form  expressions  for  the  sensitivities  of  X  ^  ,  ij ^  and  y^ 
with  respect  to  parameters  in  A,  provided  all  j  of  A  are  distinct  (Ref  2: 
163-170).  Their  results,  in  the  terminology  of  this  paper,  follow.  By 
definition  A,  A^  ,  Uj  and  v.  satisfy  the  equations 


fl  h  =  b  aj 


v  .  A  =  v  .  A  . 
-J  -J  J 


—fn  *-■£  =  —6  — m  ~  " Em  l,  m  =  1  ,  2,.  .  .  ,  NA 


j  =  1,  2,.  .  .  ,  NA 


where  5 ,,  is  the  Kronecker  delta: 
un 
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(48) 

(49) 

(50) 


i 


1  if  £  =  m 


£m 


0  if  l  f  m 


(51) 


Differentiation  of  Eq  (48)  with  respect  to  9.  and  premultiplication 


by  Vj  yields 


v.  A  /  •  <  u.  +  A.  v.  u.  =  v.  A,  u,  +  A.  v.  u. 

J  (  )  J  J  J  J(i)  J  J(i)  J  J  J  J(i) 


i  =  1 ,  2, .  .  NPA 


(52) 


where  NPA  equals  the  number  of  parameters  in  A.  Solving  Eq  (52)  for  A. 

J(i) 

yields 


A-  =  REIGV(£,j) 

°(i) 

*  EIGV(m, j ) 

(53) 

where 

REIGV  =  |  v 1  v 2  ’  ' 

^na] 

(54) 

EIGV  =  [  ^  u2  •  ' 

^NA  ] 

(55) 

and  9.  appears  in  the  l, m  location  of  A. 

Through  similar  operations  on  Eqs  (48)  and  (49)  and  by  introducing  a 


matrix  H  with  components 

h  l,rn  REIGV(l,k)  *  EIGV(m,j) 

k,j  '  Xj  -  Xk  (56) 


REIGV(l,  j )  *  EIGV(m,k) 


(57) 


the  eigenvector  and  reciprocal  eigenvector  sensitivities  may  be  written 
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h 


% 


(i) 


NA 

Z 

k=l 

Mj 


k,j 


(58) 


-j 


(i) 


NA  i,  m 

r  /-  * 

^  i  k 

k=l  J,K 

k^j 


4 


(59) 


Other  Vector  Sensitivities.  Since  it  is  assumed  that  appears 
linearly  in  B,  C  and  ,  the  partial  derivatives  of  these  vectors  with 
respect  to  contain  all  zero  elements  except  in  the  location  of  9^  , 
where  unity  appears.  The  computational  aspects  of  this  assumption  are 
discussed  in  the  next  chapter. 

Singular  Value  Decomposition  and  G.  At  this  point,  all  the  equations 
necessary  for  the  calculation  of  S  have  been  derived.  It  can  be  shown 
that  in  order  for  S  to  have  rank  NP  (  a  necessary  condition  for  a  solution 

E  FI  and  G 


of  the  identification/estimation  problem  to  exist),  both 
must  have  rank  NP  (Ref  1  1:  244;  18:  91).  The  requirement  on  the  ranks  of 
these  matrices  leads  directly  to  time,  input,  and  structural  conditions 
on  identif iabi 1 i ty .  Only  the  structural  condition  on  identifiabil ity,  i.  e., 
that  G  have  rank  NP,  is  considered  in  the  following  discussion. 

The  rank  condition  on  G  may  be  shown  to  be  satisfied  by  performing  a 
singular  value  decomposition  on  G.  As  in  Eq  (10),  the  number  of  samples 
K  equals  or  exceeds  the  total  number  of  parameters  NP.  Then  there  will  be 
NP  singular  values  if  G  has  rank  NP.  The  singular  value  decomposition  of 
G  also  provides  very  important  information  on  the  "structural"  aspects  of 
the  problem:  Calculating  the  condition  number  of  G  (See  Eq  (23).)  can 
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give  the  designer  a  good  indication  of  just  how  well  the  model  defined  by 
Eqs  (1)  through  (3)  describes  the  system.  This  capability  is  important, 
for  the  parameter  estimation  can  only  be  as  good  as  the  model.  Also,  if 
the  condition  number  of  G  indicates  that  the  model  is  good,  but  the  con¬ 
dition  number  of  S  shows  that  the  estimation  is  not,  the  designer  knows  to 
change  either  the  sample  spacing  and/or  the  input  and  number  of  samples 
to  improve  the  estimation.  A  more  detailed  description  of  this  use  of  the 
condition  number  in  the  design  phase  is  presented  in  the  next  section. 

Points  for  the  Modal  Method  in  the  Design  Phase 

The  design  phase  of  the  parameter  identification/estimetion  problem 
involves  evaluating  the  sensitivity  matrix  and,  if  necessary,  changing 
some  part  of  the  total  system  (i.  e. ,  input,  sample  spacing  and  system 
model ) . 

In  order  to  calculate  the  sensitivity  matrix  via  the  sensitivity 
system  method,  the  designer  must  form  the  ( NA  +  1)  *  (NP  +  1)  eouations 
of  the  sensitivity  system  (Eqs  (29)  and  (30))  from  the  system  equations 
(Eqs  (1)  through  (3))  each  time  a  design  change  is  made.  Although  the 
integration  is  made  simple  by  using  available  software  packages,  setting 
up  the  sensitivity  system  can  be  a  tedious  chore.  Using  the  interactive 
computer  program  which  implements  the  modal  algorithm,  however,  all  the 
designer  need  do  is  input  the  matrices  and  vectors  of  Eqs  (1)  through  (3) 
and  view  the  results  as  the  computer  prints  them  out.  If  the  results  are 
less  than  desirable,  the  matrices  and  vectors  are  easily  changed  without 
having  to  rewrite  the  program. 

In  addition,  the  sensitivity  system  method  gives  no  indication  of  how 
well  the  model  describes  the  system.  By  contrast,  in  the  modal  algorithm 
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the  time-  and  input-dependent  parts  are  separate  from  that  part  of  the 
sensitivity  matrix  which  depends  on  the  model.  See  Eq  (33).  Because  of 
this,  the  designer  can  look  at  the  condition  number  of  6  and  decide 
whether  the  model  is  adequate.  If  he  decides  it  is,  but  the  sensitivity 
matrix  is  too  ill-conditioned,  the  designer  may  conclude  that  some  aspect 
of  the  input  (sample  spacing  A  or  number  of  sample  times  K)  is  at  fault 
and  amend  the  situation.  Since  the  program  listed  in  Appendix  B  is  inter¬ 
active,  any  parts  of  the  system  are  easily  changed--wi thout  having  to  make 
up  a  new  set  of  differential  equations.  Examples  are  presented  in  Chapter  IV. 

Summary 

Both  the  sensitivity  system  method  and  the  modal  method  of  calculating 
the  sensitivity  matrix  S  of  a  linear,  time-invariant  control  system  may  be 
used  in  identifying  and  estimating  system  parameters  via  quasilinearization. 
The  theory  of  quasilinearization  and  the  two  methods  of  calculating  S  have 
been  explained,  and  some  of  the  advantages  of  the  modal  method  over  the 
sensitivity  system  method  have  been  noted.  i 

The  singular  value  decomposition  and  the  insights  it  gives  into  system 
identification  and  parameter  estimation  have  oeen  touched  upon.  In  addition, 
a  structural  condition  on  identi f iabil i ty  has  been  formed  by  combining  the 
technique  of  singular  value  decomposition  with  the  modal  method. 

From  an  analysis  viewpoint,  the  new  algorithm  seems  superior  to  the 
sensitivity  system  method.  What  remains  to  be  seen  is  how  easily  the 
algorithm  may  be  implemented  on  the  computer.  The  computational  load  is 
the  subject  of  the  next  chapter. 
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III.  Computations 

The  modal  algorithm  for  calculating  the  output  sensitivity  matrix 
was  developed  in  the  last  chapter,  and  a  listing  of  the  interactive  program 
implementing  the  algorithm  appears  in  Appendix  B.  The  interactive  program 
proves  useful  to  the  system  designer  in  both  the  experimental  design 
phase  and  in  the  estimation  phase  of  the  system  identification/estimation 
task.  Using  the  interactive  program  in  the  experimental  design  phase, 
the  designer  may  analyze  alternative  models,  alternative  sample  spacings 
and  alternative  inputs  using  experimental  data  before  the  actual  parameter 
estimation  task.  During  the  second  phase,  the  estimation  phase,  the 
computational  efficiencies  of  the  program  play  a  fundamental  role  because 
in  the  iterative  quasilinearization  technique  (explained  in  Chapter  II) 
the  sensitivity  matrix  must  be  calculated  at  each  iteration. 

Some  of  the  computational  efficiencies  of  the  modal  algorithm  are 
noted  in  the  first  section  of  this  chapter.  In  the  second  section,  the 
program  implementing  the  modal  algorithm  is  explained;  and  its  comput¬ 
ational  loading  is  tabulated.  Then  variations  of  the  basic  assumptions 
of  the  algorithm,  such  as  the  nonlinear  appearance  of  the  parameters  in 
the  system  vectors  B,  £,  and  x^ ,  and  their  effect  on  the  program  are 
discussed . 

Computational  Efficiencies  of  the  Modal  Algorithm 

In  this  section  some  of  the  computational  efficiencies  of  the  computer 
program  listed  in  Appendix  B  are  summarized. 

The  modal  algorithm  uses  exponentials,  multiplications,  additions, 
and  transfers  from  one  matrix  to  another.  The  program  takes  advantage  of 
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the  sparsity  of  the  sensitivity  matrices  and  vectors  such  as  EIGV^.j  and 
C(-).  Rather  than  multiplying  zero  as  well  as  nonzero  matrix  or  vector 
elements,  the  computer  program  calculates  only  values  affected  by  the 
parameters.  For  example,  since  P  ■  is  assumed  to  appear  linearly  in  B_,  £ 
and  x^,  the  sensitivity  vectors  Bj.^,  C ^  ^  and  x^  .  have  all  zero  elements 
except  in  the  location  of  0.,  where  unity  appears.  Therefore,  scalars 
such  as  ik  do  not  require  the  multiplication  of  two  vectors,  but 

involve  simply  a  transfer  of  the  appropriate  element  from  the  vector  u. 


to  the  corresponding  location  in  C j  ^  .  This  example  is  discussed 

further  later  in  this  chapter.  The  transfers  just  mentioned  decrease 
computation  time  and  do  not  cause  roundoff  errors.  Another  case  in  which 
advantage  is  taken  of  sparse  matrices  occurs  when  x^  is  the  zero  vector: 
Computation  time  is  decreased  by  not  multiplying  any  vectors  or  matrices 
associated  with  xfl .  Instead  of  using  Eq  (33)  to  calculate  E,  since  G,. 
is  identically  zero  when  x^  is  0  (See  Eq  (47).),  S  becomes  F  *  3ZS. 

Another  "high-quality"  feature  of  the  program  is  the  use  of  readily 
available  software  packages  to  compute  the  eigenvalues  and  \ectors  (Ref  o: 
EIGRF)  and  the  singular  values  and  vectors  (Ref  5:  LSVDF).  'he  singular 
value  decomposition  package  requires  that  the  matrix  to  oe  decomposed 


be  real  (Ref  5:  LSVDF).  However,  the  matrix  to  be  decomposed 


'zs 


,  G  or 


S)  is  expressed  in  complex  notation.  Therefore,  the  vector  of  eigenvalues 
is  checked.  This  entails  checking  NA  rather  than  2  NA  *  NP  (for  Gz<;), 

4  NA  *  NP  (for  G)  or  K  *  NP  (for  S)  values.  This  check  wil1  be  explained 
in  the  nex -  section . 


The  Program 

In  this  section  the  program  implementing  the  modal  alc'rithm  will  be 

explained.  Stress  is  placed  on  the  interaction  between  the  user  and  the 
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computer:  The  user  inputs  information;  and  the  computer  calculates  with 
the  information,  stores  it,  asks  for  more,  prints  out  results,  and  so  on. 
Once  the  program  has  been  called,  the  user  has  the  option  of  seeing  an 
explanation  of  the  matrices  and  vectors  used  in  the  algorithm.  If  the 
user  opts  not  to  see  the  explanation,  he  will  be  prompted  for  the  neces¬ 
sary  data  to  calculate  first,  the  structure-dependent  part  and  second,  the 
input-  and  time-dependent  parts  of  S. 

Structure-Dependent  Part  of  S.  The  user  inputs  the  dimension  of  A 

(‘.'A);  the  number  of  parameters  in  A  (NPA),  x^  (NPX)„  B  ( N  P  B ) ,  and  C 

(NPC);  and  the  total  number  of  system  parameters  (NP).  The  user  is  then 

asked  if  x  equals  0.  If  it  does,  the  computer  sets  x  to  0;  if  not, 

— o  -  -o  — 

the  computer  asks  for  x^.  The  system  plant  matrix,  input  and  output 
vectors  and  the  matrices  storing  the  numbers  and  locations  of  parameters 
in  A,  x  ,  B  and  C  are  requested  and  entered.  The  matrices  are  read  in 
by  columns,  as  mentioned  earlier.  The  matrices  containing  the  parameter 
numbers  and  locations  are  IA,  IB,  IC  and  IX.  An  example  of  the  contents 
of  IA  follows:  If  ?.  appears  in  the  (<  ,r)  location  of  the  A  matrix,  the 

parameter  number  i  is  stored  in  location  IA(l.LPA).  LPA=  1.2 . NPA. 

The  row  address  *  and  column  address  n  of  the  oerameter  are  stored  in 
locations  I A ( 2 . L  PA ;  and  IA(3,LPA).  respectively.  Fora  single-input 
single-output  system,  8  and  C  are  one-dimensional  arrays.  Therefore,  IB 
and  IC,  as  well  as  IX,  have  dimension  2  x  NPB,  NPC  or  NPX.  The  parameter 
numbers  are  stored  in  the  first  row.  and  the  vector  locations  of  the 
parameters  are  stored  in  the  second  row  of  IB.  IC  and  IX. 

I  ms  L  subroutine  EIGRF  (EIGCC  for  complex  A)  computes  the  eigenvalues 
and  eigenvectors  of  A.  The  eigenvalues  are  stored  in  ihe  vector  EIG,  and 
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the  eigenvectors  are  stored  by  columns  in  the  matrix  E1GV,  as  in  Eq  (55). 
The  program  then  checks  for  repeated,  zero  and  complex  eigenvalues.  If 
A  has  a  repeated  eigenvalue  or  ar  eigenvalue  equal  to  zero,  the  program, 
stops;  and  the  user  must  remodel  the  system  model  so  that  A  has  NA 
distinct,  nonzero  eigenvalues.  The  matrix  A  must  have  distinct  eigenvalues 
or  the  denominators  in  Eqs  (56)  and  (57)  will  be  zero.  The  eigenvalues 
must  be  nonzero,  or  the  denominators  of  the  constants  C4  and  C5  will  be 
zero.  See  Eqs  (42)  through  (45).  The  effect  of  the  check  for  complex 
eigenvalues  is  explained  in  the  discussion  of  IV.SL  subroutine  LSVDF  later 
in  this  section. 

To  form  the  m,atrix  of  reciprocal  eigenvectors,  I  MS  L  subroutine  LEQT1C 

is  called  to  invert  the  transpose  of  the  matrix  of  eigenvectors.  Each 

column  of  the  matrix  REIGV  is  a  reciprocal  eigenvector,  as  in  Eq  ^54). 

"-er.  v.  >  v.  8  and  C  u.,  the  vectors  which  do  not  depend  on  the  parameter 

are  formed.  The  user  may  choose  to  see  these  if  he  wishes. 

're  program  then  calls  its  own  subroutine  SENS  to  form  the  sensitiv- 

•*.  ies  cf  '•  .  .  u . .  v.,  x  ,  B  and  C  with  respect  to  each  parameter  .  Since 
j  j  j  -o  -  •  -  i 

..  appear  more  t^an  once  in  the  system  matrix  a-d  vectors,  the  $  e-  s  i  - 


t  i .  -ties  with  respect  to  are  formed  first;  then  sensitivities  with 
re:  meet  to  and  so  on  to  -,ID  are  formed.  This  involves  checking 

t-e  matrices  IA,  IB,  IC  and  IX,  which  hold  the  number  and  locations  of 
each  pere-eter,  at  nest  2  *  NP  times.  Since  the  sensitivity  vectors 

-  ,  .  .  and  have  zero  or  unity  elements,  the  values  of  u .  and  v. 

•  l  .•  i ;  ~  -i  -J 

:  c  m  >-  e  s  r  c  n  d  i  r  c  to  unity  elements  in  *ne  sc-rsitivity  vectors  are  transferred 

rr.  the  expropriate  locations  cf  C  ( ^ ,  u..  v.  x^^  and  Bj.^.  The 

qua- *  i *  •  f--  .  ,  u  and  v  .  are  cor  put  ed  via  methods  proposed  tv 

J  ( i ) 


/ 


Crossley  and  Porter  as  described  in  Chapter  II.  When  the  plant  matrix  A 

is  in  the  diagonal  canonical  form  and  the  parameters  are  the  elements 

along  the  main  diagonal  of  A,  the  eigenvalues  of  A  are  the  parameters 

themselves:  and  both  EIGV  and  RE3GV  equal  the  identity  matrix.  Therefore, 

subroutine  SENS  could  be  simplified  for  the  special  case  of  the  diagonal 

canonical  form.  But  this  is  a  topic  to  be  left  to  further  research. 

If  x  is  the  zero  vector,  G  .  is  identically  zero:  and  subroutine 
—o  z  1 

G’-iATR  is  called  to  form  G  If  x  is  not  the  zero  vector,  subroutine 

ZS  -0 

3MATR  forms  G^ ■  and  then  Gz$  .  G  is  then  formed  as  in  Eq  (47).  As  may  be 

seen  from  this  equation,  G  is  complex  if  any  eigenvalues  or  eigenvectors 

of  A  or  any  of  the  system  vectors  ,  B  or  C  are  complex.  This  is  the 

point  at  which  the  eigenvalue  check  mentioned  earlier  enters;  for  in 

order  to  perform  a  singular  value  decomposition  on  a  co~pl ex  matrix,  it 

~jst  be  put  into  a  real  number  format  (Ref  5:  LSVDF) .  .herefore,  if  the 

eigenvalue  is  real,  the  associated  row  of  the  complex  matrix  (say  Gzs)  is 

transferred  to  a  real  matrix  (say  G  ).  If  the  ith  eigenvalue  is  the 

“real 

first  or  second  of  a  complex  conjugate  pair,  the  reel  or  complex  part  of 

f.e  ith  row  of  G  is  transferred  to  the  ith  or  i  1st  row  of  G 

2S  “real 

respectively.  For  example,  if  (in  complex  notation, 

EIG  =[(1,0)  (1,1)  (1,-1)  (2,0)]  (60) 

and 

[(1,0)  (2,  0)  (3,  0)  (4,  0) 

(2,  1  )  (1,3)  (3,  -2 )  (5,-4) 

r  - 

'  (2,-D  (1,-3)  (3,2)  (5,4) 

I  (3,  0)  (2,  0)  (5,  0)  (1  .  0)  (61  ) 

U  -i 
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then 


(62) 

This  procedure  is  possible  since  I MS L  subroutine  EIGRF  beeps  complex 
conjugate  pairs  of  eigenvalues  and  eigenvectors  together  (Ref  5:  EIGRF). 

The  singular  values  and  condition  number  of  G  or  Gzs  are  then 
printed.  If  at  this  time  the  user  wishes  to  change  the  system  matrices 
and  vectors,  he  must  type  ""A"  (return)  and  begin  the  program  again.  If 

no  change  is  desired,  the  program  forms  the  time-  and  input-dependent 

parts  of  the  sensitivity  matrix. 

Time-  and  Input-Dependent  Parts  of  S.  Once  G  has  been  formed,  the 
program  prompts  the  user  to  input  the  sample  spacing  .  the  number  of 
Sc'-ples  taken  K,  and  the  discretized  input  at  each  time  t^. 

The  matrices  E  and  F  are  formed  according  to  Eas  (38)  through  (AS): 

and  if  x_^  is  not  the  zero  vector,  S  is  computed  using  Eg  (33).  If  x^  is 

the  zero  vector,  then 

S  =  F  *  G,  (63) 

L  S 

since  G  .  and.  therefore.  E  *  G  .  have  all  zero  elements.  The  user  can 
zi  zi 

perform  a  singular  value  decomposition  on  S,  if  desired.  Since  S  is  a 
real  matrix  (for  the  real  matrix  A)  but  has  complex  notation  in  the 
program.,  the  real  matrix  S  simply  contains  the  real  parts  of  the 
complex  matrix  S.  I'-'.SL  subroutine  LSVDF  nay  then  be  called  to  perform 
a  singular  value  decomposition  of  S,  and  the  condition  number  of  S  may 
oe  calculated. 


zs 


real 


12  3  4 

2  13  5 

1  3-2-4 

3  2  5  1 
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Computet ional  Load 

With  the  assumptions  oT  constant  sample  spacing  /.  and  a  piecewise 
constant  input  U,  the  computations  of  matrices  £  and  F  become  extremely 
efficient.  For  details,  see  the  development  of  E  and  F  in  Appendix  A 
'resulting  equations  are  presented  in  Chapter  II.)  and  subroutine  EFMAT 
in  the  program  listed  in  Appendix  B. 

Subroutine  GMATR,  in  which  both  G  .  and  G  are  formed,  involves 

zi  zs 

only  the  multiplications  and  additions  shown  in  the  second  half  of  Eq  (47). 
Subroutine  SENS  calculates  the  sensitivities  needed  for  forming  G,  as 
exp! a ined  earl i er  . 

The  full  computational  load  is  presented  in  Table  I.  As  shown  in 
the  teb’e,  : h e  number  of  multiplications  needed  to  compute  G  is  of  order 
’.A  *  h-A,  and  the  number  of  multiplications  needed  to  form  E  and  F  is 
or  order  *;A  *  K .  Thus  with  the  assumptions  of  constant  sample  spacing, 
oic-cewise  constant  input,  and  the  linear  appearance  of  parameters,  the 
computational  load  of  this  algorithm  is  snail. 


a r i e t i cn_s  c-f  the  Assumptions 

If  the  assumptions  on  the  sample  spacing,  input  and  parameters  are 
changed,  implementing  the  algorithm  becomes  more  difficult.  The  multi- 
input  multi-output  case  involves  more  computation  but  does  not  increase 
the  complexity  of  the  computations.  Examples  are  presented  in  the 
rcl lowing  sections. 


line  dual 


ina.  If  the  sample  spacing  is  allowed  to  vary. 


then  the  recursive  formulas  of  the  E  matrix  are  no  lonaer  valid.  For 


ex  am  pi  e  ,  e 

4- 

Vi 

e  ~  o  s  t 


does  not  necessarily  equal  e 


t 


(See  Eq  ( 38} .  )  .  and 


be  calculated  for  each  i=  1  .  2 . ND  ar.d  k  =  1 


'.ew  fcrrulas  must  also  he  derived  to  compute  the  F  matrix 
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C  ,  . 


TABLE  I 


COMPUTATIONAL  LOAD  OF  MODAL  ALGORITHM 


Code  Section  Number  of  Multiplications  Number  of  Additions 

IMSL :  EIGRF  , 


(EIGCC) 

NA 

0) 

IMSL:  LEQT1C 

NA2 

(2) 

? 

NA 

(2) 

IMSL:  LSVDF 

2K(NP)2 

(3) 

2K(NP)2 

(3) 

Vectors  independent 
of  parameters 

3NA 

3NA 

Sensitivities  of 
eigenvalues  X. 

J(i) 

Sensitivities  of 
eigenvectors  u. 

~J(i) 

and  v. 

-J(i) 

NA  *  NPA 

0 

2NA  *  NPA 

NA  *  NPA 

Sensitivities  in  x  , 

B  and  C 

0 

0 

Calculation  of 

G  .  and  G 
zi  zs 

2MA  *  (NP  +  1) 

2NA  *  NP 

Calculation  of  E  and  F 

9NA  *  { K  +  1 } 

3NA  *  (K  + 

1) 

Calculation  of  S 

G  .  f  0 
n 

4K  *  NA  *  NP 

3K  *  NP 

G  .  =  0 
zi 

2K  *  NA  *  NP 

K  *  NP 

Calculation  of  - 

1 

0 

i  i 

Notes:  (1)  Values  unavailable 

(2)  Ref  3:  1 .22 

i  (3)  Ref  3:  11.18  j 


Input  Not  Piecewise  Constant.  If  the  input  is  not  piecewise  constant, 
the  integrals  of  the  F  matrix  must  be  integrated  using  numerical  analysis 
methods.  The  equations  for  computing  F  become  more  complex,  and  their 
calculation  requires  more  computer  time. 

Nonlinear  Appearance  of  Parameters.  If  the  parameters  are  allowed 
to  appear  nonlinearly  in  the  vectors  B_,  £  and  x^ ,  the  sensitivities  of 
the  appropriate  terms  of  these  vectors  would  not  necessarily  be  unity. 
Therefore,  the  partial  derivatives  with  respect  to  the  parameters  would 
have  to  be  taken  and  multiplied  to  the  corresponding  elements  of  u_j  and  v_.. 
The  values  would  have  to  be  placed  in  the  appropriate  locations  of  £,.^  , 


v .  x  and  v  .  B , 

C  -  [  2 

f  .  v .  For  exampl  e  ,  if 

k  ‘  / 

e52  ^ 

(64) 

where  S^=  2  ,  and 

Ai  *  [  3 

4  ]' 

(65) 

then 

-(5)  =  C  ° 

2e5]  * 

’  3 

=  16 

4 

(66) 

Similarly,  if  5.  appears  nonlinearly  in  the  (•£'.  m)  location  of  A.  then 

a,  must  be  calculated  and  multiplied  to  the  right-hand  sides  of 

t,m,n 

Eas  (  53)  ,  (58)  and  (59) . 

The  Multi-Input  Multi-Output  Case.  For  a  single-input  single-output 
system  the  vectors  £  and  £  of  Eqs  (1)  and  (3)  have  dimension  NA  x  1  and 
1  x  NA,  respectively.  However,  if  NB  equals  the  number  of  inputs  and  NC 
equals  the  number  of  outputs,  B  and  £  become  matrices  B  and  C  of  dimension 
NA  x  NB  and  NC  x  NA,  respectively.  Then  the  program  LKP  would  become  a 
subroutine  to  be  called  NB  *  NC  times  from  a  new  main  program!.  For  example, 
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if  the  system  is  described  by  the  equations 


'  1 

o' 

1 

2" 

Y 

X 

- 

x  + 

0 

2 

3 

4 

U2 

. 

- 

L  c.  J 

2 

3 

*(tk} 

= 

4 

5 

then  one 

of 

t  he 

four  cal  Is 

to 

the 

(67) 


(68) 


information 


1  0 
0  2 


x  + 


(69) 


y2(tk)  =  [  4  5  ]  x(tk)  (70) 

The  "subroutine"  LKP  would  compute  four  sensitivity  natrices,  one  for 
each  input-output  pair. 


Summary 

The  use  of  the  modal  algorithm  and  the  computer  program  implementing 
the  algorithm  in  the  experimental  design  plase  and  in  the  actual  esti¬ 
mation  task  has  been  noced.  The  computational  efficiencies  of  the  program 
have  also  been  summarized,  and  the  computational  load  or  the  program  has 
been  tabulated.  In  addition,  the  effects  of  changing  some  of  the  basic 
assumptions  of  the  algorithm  have  been  mentioned  and  exemplified.  In  the 
next  chapter,  the  algorithm  will  be  verified  and  used  to  investigate 
several  experimental  design  issues. 
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IV.  Results 


Not  only  can  the  sensitivity  matrix  be  calculated  using  the  new  algo¬ 
rithm  developed  in  Chapter  II,  but  several  experimental  design  issues  can 
be  addressed.  The  design  issues  are  divided  into  two  categories:  input 
design  issues  and  structural  design  issues.  Under  the  category  of  input 
design  issues  fall  changes  in  the  sample  spacing  a,  the  number  of  sample 
times  K ,  and  the  input  U(t^).  Structural  design  issues  include  choosing 
a  zero  or  nonzero  initial  state  and  varying  the  form  of  the  mathematical 
model  by  adding  zeros  or  poles  to  the  transfer  function  or  by  choosing 
different  canonical  forms  to  represent  the  system.  Since  the  program 
implementing  the  algorithm  is  interactive,  the  user  may  easily  change  the 
input  and  structural  design  to  find  the  best  model  of  the  system  before 
the  final  experiment  must  be  performed. 

In  the  first  section  of  this  chapter,  the  algorithm  is  verified  by 
comparing  its  result  with  that  of  the  ,:starderd"  sensitivity  system  method. 
In  the  second  section,  the  experimental  design  issues  of  several  examples 
are  investigated. 

Verification  of  the  Algorithm 

To  verify  the  algorithm,  the  sensitivity  ratrices  of  two  systems  are 
formed  by  using  the  sensitivity  system  method  end  then  by  using  the  new 
al corithm. 

The  first  system  is  described  by  the  equations 


0 

1 

x  + 

A 

"1 

-2 

-2 

- 

- 

i  C  J 
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where  6  =  [0  1  10  100]  .  The  sensitivity  system  results  in  a  set  of 
six  ( N A  *  (NPA  +  NPB  +  NPC)  +1=6)  state  and  input  equations  and  five 


:  i  r 

i  y ( t  )  !l0  100  0  0  0  0 

K  I 

!y(1 j(tk)  jO  0  10  100  0  0 

!  y ( 2 ) ( t  k )  =  I  0  0  0  o  10  100 

y,3)(tk)  !  1  00  00  0 

jy(4)(tk)j  :  0  10  00  0 


Three  cases  are  considered: 


1 .  Cas  el:  x  = 
-o 


,  Ud( t k )  =  0  for  k  =  1  .  2 . 10, 


(74) 


(zero  input). 


2.  Case  II:  x  =  „  U,(t.)  =  0  for  k  =  1 ,  2,.  .  10 

—0  1  OK 

(zero  input). 

0 

3.  Case  III:  x  =  ,  U.(t.)  =  1  for  k  =  1,  2,.  .  .,  10 

“°  0  d 

(zero  state,  unit  step  input). 

The  ODE  software  package  (ReflO)  was  used  for  the  integration.  The 
sensitivity  matrices  formed  by  both  methods  may  be  seen  in  Tables  II,  III 
and  IV.  In  Cases  I  and  II,  the  input  corresponding  to  the  first  two  pa¬ 
rameters  (both  located  in  the  input  vector  B)  is  zero,  with  the  result 
that  the  first  two  columns  of  the  sensitivity  matrices  are  identically 
zero.  Therefore,  only  the  second  two  columns  appear  in  Tables  II  and  III. 

As  seen  in  Tables  II  through  IV,  the  numerical  results  of  both  methods 
are  very  close.  The  reader  might  note  that  the  new  algorithm  is  inherently 
more  accurate  than  the  sensitivity  system  calculation  since  the  errors  of 
numerical  integration  are  not  introduced  in  the  modal  sensitivity  method. 
Unfortunately,  as  the  eigenvalues  of  A  approach  each  other,  the  sensitivity 
matrix  tends  to  become  ill-conditioned  for  both  methods.  The  output 
response,  however,  may  be  relatively  insensitive  to  changes  in  the  parameters 
even  wiien  the  parallel  eigenvectors  drive  the  sensitivity  matrix  to  infinity. 
To  avoid  this  situation,  the  system  may  be  written  in  a  block  diagonal  form. 
Then  the  eigenvalues  of  A  lie  along  the  main  diagonal  of  A,  and  the  eigen¬ 
vector  and  reciprocal  eigenvector  sensitivities  are  identically  zero. 

This  not  only  simplifies  the  formation  of  the  sensitivity  matrix,  but  also 
alleviates  the  problem  of  nearly  equal  eigenvalues .  This  problem  and  its 
solution  are  discussed  in  detail  by  Reid  and  Palmer  (Ref  16:  21-33). 

Another  system  used  to  verify  the  modal  algorithm  is  described  by  the 
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tabu:  ii 


0.9  .571204708.3421  -.6369939102161  .571204708404  -.6369539102258 


TABLE  III 


.3184769549788  -.06574920184484  .3184769551129  -.06574920182172 


TABLE  IV 


0.9  -35.406806703133.9916719684 .214397645857 .31  84 76955098f  35 . 406806692  33 .991 671 9692  .214397645798  .31  84769551  1  ? 


equations 


(75) 

(76) 


When  the  nominal  values  of  the  parameters  are  f_  =  -1  -2 

the  system  becomes 
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and 
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The  zero  state,  unit  step  case  sensitivity  matrices  formed  by  the  two 
methods  appear  in  Table  V. 


Design  Issues 

To  investigate  structural  design  issues,  the  system  transfer  function 

2  2 
G(S)  = - - - —  =  -7 - 

(s  -1-  1 )  (s  +  2)  s  +  3s  +  2  (79) 

is  expressed  in  both  the  phase  variable  and  the  diagonal  canonical  forms, 
the  zero  state  part  of  the-G  matrix  (G  )  is  formed  for  both  cases,  and 
their  condition  numbers  are  calculated  and  compared.  Ey  comparing  the 
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condition  n utbe rs ,  a  conclusion  may  be  drawn  as  to  which  form,  is  better 


conditioned  and,  thus,  preferable  over  the  other  for  representing  the 
system  in  the  identification/estimation  problem.  Then  one  of  the  forms 
is  chosen  to  investigate  two  of  the  input  design  issues  mentioned  earlier. 
Ey  calculating  and  plotting  the  inverse  condition  numbers  of  successive 
trials,  the  effects  on  the  output  sensitivity  of  changing  /•  and  K  may 

^  P>  P  ^  1 6c! 

Structural  Design  Issues .  "he  o('ase  variable  form  of  Eg  (79)  is 


The  diagonal  form  is 


w^e-e  I  =  |-1  -2  2  -2|  . 

The  input  is  of  no  interest  for  the  purposes  of  this  investiaation; 

1 

therefore,  only  the  inverse  cf  the  condition  number  is  recorded.  As 

k 

shown  in  table  VI,  the  diagonal  form  has  a  better  condition  number  than 
the  chase  variable  form  has.  'his  result  is  typical,  implying 'that  the 

<1 


TABLE  VI 
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Note:  (1)  No  zeros  in  transfer  function. 


output  sensitivity  matrix  using  the  diagonal  form  remains  better  conditioned 
as  model  parameters  change  than  does  that  using  the  phase  variable  form.  The 
reason  for  its  better  conditioning  is  that  the  diagonal  form  gives  a  balance 
of  the  sensitivity  among  the  various  parameters,  yielding  a  better  cond¬ 
itioning  of  the  inherent  linear  equation  problem  of  quasi  1 inearization . 

On  the  other  hand,  the  phase  variable  form  may  give  a  high  sensitivity  to 
some  parameters  and  low  sensitivity  to  others.  This  imbalance  in  the 
sensitivity  yields  a  poorly  conditioned  problem  which  has  large  errors  in 
certain  "directions,"  as  discussed  in  Chapter  II. 

To  further  compare  the  canonical  forms  and  to  see  the  effect  of  adding 
zeros  to  the  transfer  function  of  a  system,  a  zero  is  added  at  -2.005 
(close  to  the  pole  at  -2)  and  than  at  -10.01  (far  away  from  both  poles). 

In  both  instances,  the  condition  number  of  the  diagonal  form  is  better 
than  that  of  the  phase  variable  form,  indicating  that  the  diagonal  form 
is  superior  in  estimating  parameters.  With  the  addition  of  zeros,  the 
condition  of  both  forms  becomes  worse  (See  Table  VI.),  indicating  that  the 
addition  of  zeros  to  the  transfer  function  makes  parameter  estimation  more 
difficult.  Another  revelation  of  Table  VI  is  that  placing  a  zero  over  a 
pole  may  lessen  the  accuracy  of  the  parameter  estimation  enough  to  make 
it  totally  unreliable. 


Input  Design  Issues .  The  investigation  of  input  design  issues  includes 
observing  how  chancing  sample  spacing  L  ,  number  of  samples  taken  K,  or  the 
input  U ^ (ty)  affects  the  output.  In  this  section  the  phase  variable  form 
of  Eq  (79)  is  chosen  since  its  condition  number  varies  mere  than  that  of 
the  diagonal  form,  as  nay  be  sc  ole  VI.  Therefore,  any.  correlation 

between  the  condition  number  cnd  changes  in  L  or  K  should  be  more  apparent 
with  this  form. 
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1 

Figure  1  shows  how  varies  with  when  10  samples  are  taken  for 
the  zero  state,  unit  step  input  case.  With  K  held  constant,  the  condition 

x 

number  improves  (  approaches  1)  as  the  sample  spacing  is  increased,  and 
then  deteriorates  when  the  sample  spacing  becomes  too  large.  The  first 
trend  may  be  surprising,  but  it  only  means  that  the  whole  interval  K  *  L 
over  which  the  sample  outputs  are  taken  becomes  larger;  and  more  infor¬ 
mation  is  gained.  As  the  interval  becomes  too  large,  however,  not  enough 

information  is  taken  between  sample  times  to  learn  the  true  output  response 

1 

This  is  reflected  by  the  decrease  in  —  . 

K 

Out  of  Figure  1  comes  the  unanticipated  discovery  that  for  each  number 
cf  sample  times,  there  is  an  optimum  sample  spacing  (and,  therefore,  an 
optimum  interval)  over  which  to  observe  the  output  response.  In  the  follow 
inc  figures,  other  systems  are  examined  to  find  the  optimum  sample  spacing 
fcr  a  given  number  of  samples. 

Figure  2  shows  the  results  for  the  system 

0  1  1  F°" 

x  +  U 

r.  2  i 

v  *1  ~  9  1  (  Q,£  ) 

y ( t k )  =  e3  e4  x(tk)  (85) 

where  6  =  -0.01  -2  0.01  0  ]  and  where  10  samples  are  taken  with 

a  unit  step  input.  The  results  for  the  system  of  Eos  (80)  and  (01)  with 
!  =  sin  t  is  shown  in  Figure  3. 

Figure  4  shows  the  relationship  of  the  inverse  condition  number  with 
sa~:le  spacing  for  the  system 
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0  1 


0  0 


0 


I 


I 


u 

(86) 

yUk)  f  e5  e6  07  08  ] 

where  6=  [  -5.05  -7.06  -7.21  -6.2  0.055  1.15  1.00  0.00  ]  "  with 

a  unit  step  input,  zero  initial  state  and  10  samples  taken.  The  curve  for 
the  fourth  order  system  is  steeper  around  the  optimum  point  than  previous 
curves.  Therefore,  it  seems  that  determining  the  optimum  sample  spacing 
during  the  design  phase  is  more  important  as  the  system  becomes  larger 
(i.  e. ,  as  the  number  of  parameters  increases). 

The  relationship  of  the  inverse  condition  number  with  the  number  of 
samples  taken  when  the  sample  spacing  is  held  constant  for  the  system  of 
Eqs  (80)  and  (81)  is  shown  in  Figure  5.  A  unit  step  input  is  used.  As 
expected,  for  constant  A  =  0.1  ,  an  increase  in  the  number  of  samples 
taken  results  in  a  better  condition  number  of  the  sensitivity  matrix  until 
a  point  is  reached  after  which  taking  more  samples  does  little  good. 

The  results  for  the  same  system  where  U  =  sin  t  may  be  seen  in 
Figure  6.  Figure  7  shows  the  results  for  the  system  of  Eqs  (84)  and  (85) 
with  a  unit  step  input  and  A  =  0.1  . 

In  Figure  8,  the  relationship  between  the  inverse  condition  number  of 
S  and  the  number  of  samples  taken  for  three  different  sample  spacings  is 
shown.  Using  this  type  of  graph,  the  designer  may  choose  an  optimum 
combination  of  sample  spacing  and  number  of  samples.  As  in  Figure  8,  the 
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4 


design  engineer  may  learn  that  the  optimum  number  of  samples  is  quite  large. 
For  A  =0.1  ,  for  example,  the  optimum  number  of  samples  (greater  than  35) 
overflows  the  field  length  of  the  interactive  terminal.  While  the  program 
dimensions  may  be  increased  and  the  program  run  on  a  batch  job,  this  defeats 
the  purpose  of  having  an  interactive  program. 

Summary 

In  the  first  section  of  this  chapter,  the  modal  method  was  verified  by 
comparing  its  results  with  those  of  the  "standard"  sensitivity  system 
method.  The  modal  method  has  been  used  to  compare  the  diagonal  and  phase 
variable  canonical  forms  and  to  get  an  idea  of  the  effect  on  the  output 
sensitivity  of  adding  zeros  to  the  system  transfer  function.  In  addition, 
the  condition  number  of  the  sensitivity  matrix  has  been  shown  to  improve 
by  increasing  the  number  of  samples  taken  up  to  the  point  where  taking 
more  samples  has  little  effect  on  the  condition  number.  And  it  has  been 
shown  that  by  varying  the  sample  spacing  in  trial  runs  of  the  program,  an 
optimal  sample  spacing  may  be  found  for  a  specified  number  of  samples. 
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V.  Conclusions  and  Recommendations 


for  Further  Research 


The  modal  method  derived  in  Chapter  II  has  been  verified  by  compering 
its  results  with  those  of  the  standard  "sensitivity  system"  method.  In 
Chapter  III,  the  computational  load  of  the  software  implementing  the  modal 
algorithm  was  analyzed.  Structural  and  input  design  issues  have  teen 
investigated  in  Chapter  IV. 

The  following  conclusions  may  now  be  made: 

1.  The  modal  method  is  computationally  more  efficient 
than  the  standard  "sensitivity  system"  method.  The 
formation  of  G  requires  on  the  order  of  NA.  *  NPA 
multiplications,  while  the  E  and  F  matrices  require 
on  the  order  of  NA  *  K  multiplications. 

2.  The  modal  method  is  inherently  more  accurate  than 
the  sensitivity  system  method  since  it  eliminates 
the  numerical  proble":  of  integrating  "stirf"  dif¬ 
ferential  equations. 

3.  The  use  of  the  singular  value  decomposition  in  con¬ 
junction  with  the  new  algorithm  leads  to  a  structural 
condition  on  parameter  identi f iabi 1 i ty .  "his  is  a 
significant  contribution  of  the  algorithm  since  the 
structural  condition  aids  in  the  investigation  of 
canonical  forms  and  in  the  overall  evaluation  of 
system  models. 

4 .  For  all  cases  in  Chapter  IV,  the  sensitivity  matrix 
of  the  diagonal  form,  is  better  conditioned  than  is 
that  of  the  phase  variable  canonical  form.  This 
implies  that  the  diagonal  canonical  form  is  more 
accurate  in  estimating  system  parameters  than  is 
the  phase  variable  form. 

5.  The  interactive  program  implementing  the  modal  method 
is  very  helpful  to  the  designer  in  both  the  exper¬ 
imental  design  phase  end  the  actual  parameter  esti¬ 
mation  task.  The  program  may  easily  be  extended 
into  the  iterative  algorithm  described  in  Chapter 

II  to  perform  the  entire  parameter  estimation  task. 


6.  For  a  constant  nurber  of  sample  times,  an  optimum 
sample  spacing  may  be  found. 

7.  for  a  constant  sample  spacing,  taking  more  samples 
increases  the  accuracy  of  the  parameter  estimation, 
until  a  point  is  reached  after  which  taking  more 
samples  has  little  effect  on  the  accuracy. 

Although  the  following  areas  have  not  been  researched  fully,  the 

available  software  allows  their  investigation: 

1.  Input  design  optimization. 

Optimum  sample  spacing  and  number  of  samples,  both 
separately  and  in  conjunction  may  be  determined. 

In  addition,  discrete  inputs  other  than  the  unit  step 
may  be  used  (Ref  S). 

2.  Initial  state  optimization. 

States  other  than  the  zero  initial  state  may  lead  to 
better  accuracy  of  the  parameter  estimation. 

3.  Nth  order  reduction. 

The  relationship  between  the  estimation  accuracy  of 
the  identifiable  parameters  and  state  model  order 
reduction  techniques  ~ay  be  investigated  (Ref  £). 

4.  Canonical  forms  anc  parameters. 

Canonical  forms  anc  syste-  parameters  other  than  the 
diagonal  and  phase  variable  forms  investigated  im 
this  paper  may  be  chosen  to  rebel  the  system.  ~r<e 
special  block  diagonal  form  discusse  in  Crapters 
III  and  IV  is  one  example  (Ref  17). 

5.  Numerical  problems  of  close  eigenvalues. 

The  growing  output  sensitivities  occurring  when  she 
eigenvalues  of  the  ulant  matrix  A  approach  each  other 
(See  Ecs  (58)  and  (55).)  should  be  given  close  atten¬ 
tion.  This  matter  is  discussed  more  fully  by  Reid 
and  Palmer  (Ref  17:  21-33). 


56 


Bib! iography 


1.  Astrom,  K.  J.  and  P.  Eykhoff.  "System  Identification--A  Survey," 

Automatica ,  7_:  123-162  (1971  ). 

2.  Crossley,  T.  R.  and  B.  Porter.  "Eigenvalue  and  Eigenvector  Sensi¬ 
tivities  in  Linear  Systems  Theory,"  Int.  J.  Control  ,  1_0  (2):  163- 

170  (February  1969). 

3.  Dongarra,  0.  J.  et  aj_.  LINPACK  Users '  Guide.  Philadelphia:  Society 
for  Industrial  and  Applied  Mathematics,  1579. 

4.  Gupta,  N.  K.  and  R.  K.  Mehra.  "Computational  Aspects  of  Maximum 

Likelihood  Estimation  and  Reduction  in  Sensitivity  Function  Calcu¬ 
lations,"  IEEE  Transactions  on  Automatic  Control,  1 9  (6):  774-783 

(December  1 974) . 

5.  IMSL  Reference  Manual .  Houston,  Texas:  IKSL,  Inc.,  1979. 

6.  Kumar,  K.  S.  and  R.  Sridhar.  "On  the  Identification  of  Control 

Systems  by  the  Quasi-Linearization  Method,"  IEEE  Transactions  on 
Automatic  Control,  22:  242-246  (April  1977). 

7.  Lav, 'son,  C.  L.  and  R.  J.  Hansen.  Sol vi no  Least  Souare  Problems. 
Englewood  Cliffs,  New  Jersey:  Prentice-Hall,  1974. 

c.  McClendon,  J.  R.  Model  Order  Reduction  Using  the  Balanced  State 
Representation :  Theory ,  Applications  and~  Interactive  Software 
Inpl  ementation .  MS  Thesis.  W  r  i  g  h  t  -  Pa  tte  r  son  A  F  B ,  Ohio":  Air  Force 
Institute  of  Technology,  December  1979. 

9.  Mehra,  R.  K.  "Optimal  Input  Signals  for  Parameter  Estimation  in 
Dynamic  Systems--Survey  and  New  Results,"  IEEE  Transactions  on 
Automatic  Control  ,  1 9  (6):  753-768  (December  1974)". 


10.  Nikolai,  Paul  J.  and  Donald  S.  Clemm.  Sol ution  of  Ordinary  Pi fferential 
Equations  on  t'ne  CPC  6600/CYBER  74  Processors .  Kright-Patterson  AFB, 
Ohio:  Air  Force  Flight  Dynamics  Laboratory,  1977. 

11.  Reid,  J.  G. ,  et  al .  "Algebraic  Representation  of  Parameter  Sensitiv¬ 
ities  in  Linear,  Time-Invariant  Systems,"  Journal  of  the  Frankl i n 

Insti tute,  301  (1  and  2):  123-141  (January  and  February  T976) . 

12.  Reid,  J.  G.  "Structural  Ider.tifiability  in  Linear  Time- Invariant 

Systems,"  IEEE  Transactions  on  Automatic  Control.  22:  242-2^6 

(April  1  977TT~ 

13.  Reid,  J.  G.  Lecture  materials  destributed  in  EE  5.10,  Linear  Systems 

Analysis  and  Digital  Computation  Methods.  School  of  Engineering,  Air 

Force  Institute  of  Technology,  Wright-Pattersor,  A^B,  Ohio,  1975. 


14. 


Reid,  J .  G .  A  New  Paral  lei  Identi fication  A1 gori thm  for  Linear  T  i  me  - 
Invariant  Systems:  Preliminary  Results.  School  of  Engineering,  Air 
Force  Institute  of  Technology,  Wright-Patterson  AFB,  Ohio,  1977. 

15.  Reid,  J.  G.  " Ident i f iabi 1 i ty  and  Experimental  Design  Issues  in  Linear 

System  Identification,"  ASME  Journal :  1  -  10  (December  1  979), 

16.  Reid,  J.  G.  School  of  Engineering,  Air  Force  Institute  of  Technology, 
(personal  correspondence).  Wright-Patterson  AFB,  Ohio,  May,  1930. 

17.  Reid,  J.  G.  and  L.  K.  Palmer.  A  New,  Highly  Efficient,  Parallel  Com- 
putation  Algorithm  for  the  Output  Parameter  Sensitivities  in  Linear 
Time-Invariant  Systems.  School  of  Engineering,  Air  Force  Institute  of 
Technology,  Wrignt-Patterson  AFB,  Ohio,  1980. 

18.  Stewart,  G.  W.  "On  the  Perturbation  of  Pseudo-Inverses,  Projections, 

and  Linear  Least  Squares  Problems,"  SIAM  Review,  1 9  (4):  634-662 

(October  1977). 

19.  Strang,  Gilbert.  Linear  Algebra  and  its  Applications.  New  York, 

San  Francisco,  and  London:  Academic  Press,  19/6. 


58 


Appendix  A 

Derivation  of  Recursive  Formulas  for  Matrices  E  and  F 


The  equations  for  calculating  the  time-  and  input-dependent  matrices 
are  derived  in  the  following  manner.  Equation  (34)  is  repeated  here: 


h  = 


r  x] tk  A2tk  XNAtk  A1tk 

e  ,  e  ,  e  ,  tk  e 


lk  e 


>'2tk 


lk  e 


XNAtk 


1  x  2NA 


(34) 


When  the  sample  times  are  all  spaced  L  apart,  t,=  L  ,  t^=  2h  ,  and  so  on 
until  tf,A=  NA  *  a  .  Then  the  first  NA  columns  of  E  may  be  computed  by 


the  recursive  formula 


El(k,f.)  =  e^*El(k-l,£) 

wi  th 

A  A 

E 1  ( 1  ,  c )  =  e  1 


(38) 


(40) 


Similarly,  the  second  NA  columns  of  E  may  be  formed  by  setting 

E3(k,£)  =  k  *  A  *  El(k,€)  (39) 

wi  th 


E3(l,£)  =  A  *  El (1 ,t) 


(41) 


In  matrix  F  the  kth  row  is  given  by 
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a  ,.a 

e  1 _ -  1 

where  C4  =  ^  .  For  k  =  2,  3,.  .  . ,  NA 

"L 


A  «kA  rku  -A/,T 

FI (k,.f)  =  e  1  I  e  U(t)  dr 
'o 


ArA  /  A»(k-1 )A  r{ k-1 ) A  -X.t 


'r  I  r 

e  e 


/\  1  /  *-  f  1 

e  U(t)  dr 


X„ki  kr  -A»t 
+  e  1  I  e  U(x)  dx 

J(k-1)A 


a A„k_, 

e  1  *  FI (k-1 ,f)  +  e  1 


~  A  *T ' 

.  L 


kA 


t  J  (k-D:. 


*  Ud(k) 


e  r  *  FI  (k-1  ,t)  +  C4  *  U  ,(k) 


(43) 


As  for  the  second  f.A  colunns  of  Eq  (35),  when  k  =  l 


F3(1 ,  C)  =  e 


¥- 


r: 

J  (A  -  t)  e  U(tj  dT 
0 

e"1'  -  1 


}'Z 

r  -  \ 


ud(l)  -  e 


¥ 


/ 


A  —  A  »x 

e  1  dx|  *  U  ,(1) 


■£  *  Ud(l)  +  e 


/, 

(-  * 


(•fT  +  De 


-  A  ,  T 

L 


*  ud<'> 
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(45) 


where  C5  = 

For  k  = 

F3(k,£) 


With  a  litt" 
te  written 


C5  *  Ud(l) 


>  A 

1  *e  * 

;  2 


X  A 


2,  3,.  .  NA,  the  second  NA  columns  of  Eq  (35)  are  given  by 


=  e 


>.„kt  kA  -Xi 

1  /  (kA  -  i)  e  1  U(t )  dT 

Jo 


ApA  /  A p ( k - 1  )A  r ( k-1  )A  ->.£t 


e  £  e  £ 


/ 


e  1  U(t)  dx 


(89) 


algebraic  manipulation,  the  first  term  of  Eq  (89)  may 


'T 


A, (k-1 )A  ,(k-l)A  \ 

;  1  J  (kA  -  (k-1 )a  +  (k-1 )A  -  i) 


-At 

e  £  U't )  di 


A  r  A  /  A  „  ( k  - 1 )  A  r  ( k  - 1 )  A 


e  f  F3(k-1 ,1)  +  A  *  FI (k-1 ,i) 


By  making  the  change  of  variable  r,  -  i  -  (k-1 ) A  ,  the  second  term  of 
Eq  (89)  becomes 


e  f  j  (kA  -  ;  -  (k-1 )a) 
o 


->.»(k-1  )a 

e  4  *  e  *  U( c- ( k-1 )a)  d: 


>-t A  ,  A 


C)  e  ^  dc  *  U,(k) 


1  r  ^ 

/ 


C5  *  ud(k) 


where  C5  is  defined  after  Eq  (45).  Combining  the  two  terms  of  Eq  (89)  gives 


F3(k,C)  =  e'r*  (F3(k-1,0  +  i  *  F1(k-1,0)  +  C5  *  Ud(k)  ( 

(Ref  17). 

The  resulting  subroutine  EFI’AT  may  be  seen  at  the  end  of  the  computer 
program  listing  in  Appendix  B. 
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Appendix  B_ 

Computer  Program  Implementing  Modal  Method  Cal cul at i on 
of  the  Sensitivity  Matrix 

The  following  is  the  listing  of  the  software  used  to  calculate  the 
sensitivity  matrix  via  the  modal  method.  The  software  is  used  on  the 
CDC  6600/CYBER  74.  Subroutines  EIGRF  (or  EIGCC),  LSVDF  and  LEQT1C  from 
the  IMSL  software  package  are  utilized  in  the  calculations  (Ref  5). 
Appendix  C  contains  the  user's  guide  to  this  program. 
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Appendix  C 

User's  Gui de  to  Interactive  Program 
Imp! ement ing  r-'pdal  Method  Calcul  at  ion  of 
the  Sensitivity  Matri x 

Program  LKP  (or  LKPC)  computes  the  sensitivity  matrix 


yO)(ti> 

y(?)(t1)  •  • 

■  y(,'.’P)^V 

y(l)(t2') 

y{2)['X2)  ■  ■ 

•  WV  i 

1 

V(i)(tK) 

y(?)(tK)  .  . 

'  y(NP)'V 

K  x  NP 

(s 

of  a  linear,  single-input  single-output,  tine-invariant  control  system  by 
decomposing  the  matrix  into  three  parts.  The  first  two  parts  depend  on 
input  and  time;  the  third  depends  solely  cn  the  structure  of  the  system 
model.  The  theory  and  algorithm  encompassed  in  the  program  are  developed 
and  explained  in  Chapters  II  and  III. 

Basic  ass'J-pticns  of  thc  algorithm  include: 

1.  The  model  is  a  linear,  tire-invariant. .  sirgl  e-incut 
sing! e-outout  system. 

2 .  The  matrix  A  and  the  vectors  3,  C.  x(t '  .  and  r  arc- 
real  (complex  for  Program  L'-'PC). 

3.  florins!  values  of  A,  B.  C,  y.  and  -  a’-"  good 
approximations  of  their  true  values. 

4.  All  are  distinct  and  nonzero. 

i 

5.  Each  •  appears  linearly  in  A.  ? .  C.  and/cr  x  ■ 

f .  •  .  rev  appear  mere  ‘her  once  and  ir  any  1cca*.''on  oJ‘ 

A .  B  .  C .  and/ o  r  x 
■  o 

?.  ?3~;le  spacing  '  is  corstan"’. 

h.  The  incut  U(t)  is  piecewise  c-mstart  . 


AO-AlOO  fi  16  AIR  FORCE  INST  OF  TECH  MRIAhT-RATTERSON  AFB  OH  SCHOO— ETC  F/#  12/1 

AN  INTERACTIVE  PROARAM  FOR  THE  CALCULATION  AND  ANALYSIS  OF  TmE  —ETC  m 
MAR  fil  L  K  PALMER 

UNCLASSIFIED  AFIT/OA/CE/81M-1  NL 


9.  If  U(t)  is  a  continuous  function,  the  program 
user  has  discretized  U(t)  over  each  inverval 
Ctk_i ,  tk]  for  k=1  ,2,.  .  . ,  K. 

10.  The  discretized  input  U{ t ^ ) ,  k=l,  2,  .  .  .  ,  K, 

is  known  exactly.  ' 

11.  The  number  of  sample  times  K  is  at  least  equal 
to  (>_)  the  total  number  of  parameters  NP. 

12.  Maximum  values  which  may  be  used  in  Program  LKP 


are : 

a.  NA  (dimension  of  A)  =  10 

b.  K  (number  of  samples)  =  35 

c.  NP  (total  number  of  parameters)  =  10 

d.  NPA,  NPB ,  NPC ,  NPX  (number  of 
parameters  in  each  of  A,  B,  C, 

and  x  )  =10 


These  values,  however,  may  be  increased 
simply  by  changing  the  dimensions  of  the 
appropriate  matrices  and  vectors  in  the 
program. 

In  explaining  the  use  of  the  program,  each  prompt  and  then  the 
possible  user  responses  are  given  in  the  order  in  which  they  occur  in 
a  computer  run.  The  computer  output  will  be  typed  in  capital  letters 
in  this  user's  guide;  the  user  responses  will  be  typed  in  lower  case 
letters.  If  at  any  time  the  user  wishes  to  terminate  a  run,  he  may 
enter  %f\  in  place  of  any  input.  (Due  to  quirks  of  the  computer,  he 
may  have  to  do  this  twice.) 

After  calling  Program  LKP,  the  following  message  will  appear: 

TITLE:  NEW  ALGORITHM  FOR  CALCULATING  THE 
SENSITIVITY  MATRIX  OF  A  LINEAR, 

SINGLE-INPUT  SINGLE-OUTPUT  (SISO) 

TIME-INVARIANT  SYSTEM 

D( X(TIME ) )/DT=A( P)*X(TIME . P)  +  B(P)*U(TIME) 

Y( DISCRETE  TIME)=C(P)*  X0( DISCRETE  TIME) 


WHERE  X0=  X(P ,T=0) 


AND  P=THETA,  VECTOR  OF  PARAMETERS  OF  INTEREST. 
A,  XO,  B,  AND  C  ARE  REAL. 


IF  A,  XO,  B,  AND  C  ARE  COMPLEX,  TYPE 
"A  (RETURN)  AND  CALL  PROGRAM  LKPC. 

DO  YOU  WANT  TO  KNOW 
HOW  TO  USE  THIS  PROGRAM? 

TYPE  I  FOR  YES,  2  FOR  NO. 

If  the  user  types  1,  definitions  of  all  the  vectors,  matrices,  and  other 
quantities  used  in  the  program  and  an  explanation  on  how  to  use  them  will 
be  given. 

After  the  example,  or  if  the  user  types  2,  the  prompt 
ENTER  NA,  NPA,  NPX,  NPB ,  NPC,  NP 

will  appear.  These  values  must  be  input  as  integers  separated  by  commas 

since  the  READ  statements  are  unformatted.  The  next  prompt  will  be 

IS  THE  NOMINAL  STATE  VECTOR 
XO  EQUAL  TO  0? 

1-YES,  2=NO . 

If  the  user  responds  with  1,  the  program  will  set  all  values  of  the 

elements  of  to  zero  and  skip  to  the  next  prompt,  if  the  response  is 

2,  the  user  will  be  asked  to  input  x  : 

—o 

ENTER  XO 

In  Program  LKP,  x  is  real  ;  in  Program  LKPC,  x  must  be  entered  in 
complex  notation. 

After  has  been  entered,  the  matrices  and  vectors  A,  IA,  B,  IB,  IX,  C, 
and  IC  will  be  requested: 

ENTER  A,  IA,  IX,  B,  IB,  C,  IC 

A,  8,  and  C  are  real  in  Program  LKP  (complex  in  Program  LKPC);  and  A  must 
be  entered  by  columns.  B  and  _C  are  one-dimensional  arrays.  IA  has  dimen¬ 
sion  3  x  NPA,  where  the  first  row  contains  the  numbers  of  the  parameters 

of  A  in  increasing  order;  the  second  row  contains  the  row  address  of  each 
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parameter,  and  the  column  address  of  each  parameter  is  stored  in  the 
third  row.  IX,  IB,  and  IC  have  dimension  2  x  NPX  (or  NPB  or  NPC) ,  where 
the  parameter  number  is  stored  in  the  first  row  in  increasing  order  and 
the  parameter  location  in  the  second  row.  The  matrices  IA,  IX,  IB,  and 
IC  are  integer  arrays  and  must  be  input  by  columns. 

Once  the  system  matrices  and  vectors  have  been  entered,  the  computer 
will  echo-print  them  along  with  the  eigenvalues  EIG,  eigenvectors  EIGV, 
and  reciprocal  eigenvectors  REIGV.  The  matrices  EIG  and  EIGV  are  calcu¬ 
lated  using  IMSL  subroutine  EIGRF  for  real  A  or  EIGCC  for  complex  A  (Ref 
6:  EIGRF,  EIGCC).  REIGV  is  obtained  by  inverting  the  transpose  of  EIGV 
via  IMSL  subroutine  LEQT1C  (Ref  6:  LEQTIC).  EIG  has  dimension  NA.  Both 
EIGV  and  REIGV  have  dimension  NA  x  NA;  and  each  column  of  EIGV  (REIGV)  is 
an  eigenvector  (reciprocal  eigenvector).  The  eigenvalues  are  then  checked. 
If  any  are  repeated,  the  message 

A  MATRIX  HAS  REPEATED  EIGENVALUES. 

PLEASE  AMEND  STATE  EQ'NS  SO  THAT 
A  HAS  DISTINCT  EIGENVALUES. 

and  the  program  will  terminate  with  the  statement 

STOP  REPEATED  E'VAL 

If  any  of  the  eigenvalues  equal  zero,  the  message 

E'VAL  CAN'T  BE  ZERO. 

will  appear,  and  the  program  will  end. 

When  all  the  eigenvalues  of  A  are  distinct  and  nonzero,  the  prompt 

WOULD  YOU  LIKE  TO  SEE  THE  CONSTANT 
VECTORS?  1 =YES  ,  2=N0. 

will  be  given.  If  1  is  entered,  v'x  ,  Cu.,  and  vtB  will  be  printed  for 
j=l  ,  2,  .  .  .,  NA.  By  typing  a  1  in  response  to 
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WOULD  YOU  LIKE  TO  SEE  THE 
PARAMETER-SENSITIVE  MATRICES? 

I =YES ,  2=N0. 

the  user  may  see  ,  Cu .  ,  vC  x  ,  v?  B,  vCx  ,  vCB,.,,  and 

J(i)  J(i )  ~J(i)~ 

£(-j)kj  j=l  >  2,  .  .  . ,  NA  and  i=l,  2,  .  .  . ,  NP.  These  values  are 
placed  in  matrices  of  dimension  NA  x  NP,  in  which  the  ith  column  contains 
the  values  associated  with  the  ith  parameter. 

If  x^  is  not  the  zero  vector,  the  structure-dependent  part  of  the 
sensitivity  matrix 


is  calculated  and  printed.  If  x  is  the  zero  vector,  however,  G  .  is 

—o  zi 

identically  zero;  and  only  Gzs  is  calculated  and  printed.  By  responding 
to  the  prompt 

DO  YOU  WANT  A  SINGULAR  VALUE 
DECOMPOSITION  OF  G=GZI  ,GZS? 

1 =  Y ES  ,  2=N0. 

the  user  may  choose  to  perform  a  singular  value  decomposition  of  G  in 

Eq  (19)  or  of  G  in  the  case  where  x^  equals  the  zero  vector.  The 
^  z  s  —o 

inverse  condition  number  —  of  G  (or  G  )  is  also  calculated.  The  user 

k  z  s 

may  then  opt  to  see  the  right  and  left  singular  vectors  or  only  the 

singular  values  and  the  inverse  condition  number  by  answering  the  prompt 

AFTER  TAKING  SING.  VAL.,  GTOT= 

UT  *  S  *  V 

TYPE  1  IF  YOU  WANT  TO  SEE 
JUST  THE  SINGULAR  VALUES  (S) 

AND  1 /CONDITION  NUMBER  OF  G  (CNG) 

TYPE  2  IF  YOU  ALSO  WANT  TO  SEE  THE 
SINGULAR  VECTORS.  ENTER  1  OR  2. 


After  the  structure-dependent  part  of  the  sensitivity  matrix  has 
been  calculated,  the  user  receives  a  message  to  enter  the  sample  spacing, 
the  number  of  discrete  times,  and  the  discretized  input: 

D=SAMPLE  SPACING.  ENTER  D. 

KIN=NUMBER  OF  DISCRETE  INPUT/OUTPUT  TIMES. 

ENTER  KIN,  UP  TO  25 

U(K),K=1 .KIN,  ARE  DISCRETE  INPUTS. 

ENTER  INPUTS. 

If  the  user  wishes  to  change  the  system  matrices  and  vectors  to 
yield  a  better  matrix  G,  he  may  enter  "%A"  in  response  to  the  preceding 
prompt  and  begin  the  program  again.  If  the  user  enters  the  values 
requested,  the  matrices  E,  F,  and  S  are  calculated.  The  E  and  F  matrices 
may  be  seen  by  entering  1  in  response  to 

DO  YOU  WANT  TO  SEE  E  AND  F? 

1 =YES ,  2=N0. 

The  sensitivity  matrix  S  is  always  printed,  and  by  typing  1  after  the 
question 

DO  YOU  WANT  A  SINGULAR  VALUE 

DECOMPOSITION  OF  SMAT? 

1 =YES ,  2=N0. 

the  singular  value  decomposition  and  inverse  condition  number  of  S  will 
be  formed.  As  before,  the  user  may  choose  to  see  the  left  and  right 
singular  vectors  along  with  the  singular  values  and  inverse  condition 
number  of  S . 

At  this  point,  if  the  user  feels  that  the  sensitivity  matrix  is  too 
ill-conditioned  or  otherwise  not  good  enough  for  his  purposes,  but  the  G 
matrix  is  satisfactory,  he  may  conveniently  change  the  input,  sample 
spacing,  and/or  the  number  of  samples  by  entering  1  in  response  to 

DO  YOU  WANT  TO  CHANGE  JUST  THE  INPUT? 

1 =YES  ,  2=N0 . 
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This  feature  of  the  program  is  particularly  helpful  when  the  user  has 
taken  too  large  a  sample  spacing  or  not  enough  samples  to  adequately 
discretize  an  input  such  as  a  sinusoid.  Then  A  and  K  may  be  changed 
without  the  necessity  of  re-entering  the  system  matrices  and  vectors. 
The  message  "D  =  SAMPLE  SPACING....11  reappears. 

Of  course,  if  the  user  answers  the  last  question  with  a  2,  the 
program  will  end  with  the  message 

STOP  NO  MORE  INPUT,  YOU  SAID. 

An  example  follows: 


TITLE:  NFI.I  ALGORITHM  FOR  CALCULATING  THE 

SENSITIVITY  MATRIX  OF  A  LI  NEAP » 

SINGLE- INPUT  SINGLE-OUTPUT  <'SISO> 

TIME- 1 N V API  ANT  *  V DTE M 

D  '■  T I  Mr')  1  DT  =  A  •  o  ,  .  T I  HE «  P  1  +  P  1  P)  ♦'_» 1  T I  ME) 

V  'III  SC  PETE  T I  ME  :■  =  •*-  ..?•)♦  XO <DI SCPETE  T I  ME) 

i.ihepf  ■/o=  x  <:p«  t*iy< 

and  P  =  THETA.  vector  of  FA  - h-'-ETEPS  OF  3NTERES 
A«  X'O,  p«  AND  C  A:E  REAL. 

7  S’  A,  :'i •  J: «  HAP  C  APE  L'DhPlE’::!«  i  YPE 
•.A  i  E T ;  l-’N  i  AND  CALL  PROGRAM  L:  P. 

DO  '.'ANT  TO  >  NDl.i 

AS1*!  TO  '.l S  E  THI I  PFD'-F^N  • 

TYC'E  1  FDP  YE* « c'  FOP  NO. 

4 7 0 0 i'i }■  CM  STORAGE  U’EP 

7  0!:  CP  SECOND":  COMPILATION  TIME 


entep  na,  npa«  up:  npi.«  npc  •  ns- 

*  M  •  1.1  (  c'  •  ^ 

i:  THE  NOMINAL  STATE  VECTOR 
■fi  EC' UAL  TO  O' 

] =YE " < £=Nn. 

1 

ENTER  A  •  I A .  I :  ‘  •  D «  I  F: «  C  «  I 

■  .  _ .  03  .  l .  -  E .  1  >  S'  •  1  •  P  •  P  <  S'  >  0  •  1  -  .  "  1  ■  ’  ■  1 '  -F  *  — 


89 


4 


AFTER  T RK I MG  X I HG . VhL . - GZS= 
UTS  ♦  S  ♦  VS 

TYPE  1  IF  YOU  1. 1  ft  NT  TD  SEE 
JUST  THE  SINGULAR  VALUES  <:'S'-» 
AND  1 -•-CONDITION  NUMBER  OF  GZS 
C'N '? Z S >  «  T Y PE  2  IF  Y Oil  A L S 0 
UANT  TD  SEE  THE  SING.  VECTORS 
ENTER  1  OR  2 
1 


1.  1  48:2066 004  .  4376285279265  ■  00577  027727 
CNGZS=  .  0  0 19151  0  0628 074 
D= SAMPLE  SPACING.  ENTER  D. 

.  1 

CI.N=NIJMEER  OF  DISCRETE  I NF'UT.-OUTPUT  TIME 
ENTER  KIN,  IIP  TO  80 
1  0 

ij  oo  .  K=  1  .  K I N  ARE  D I SCRETE  I NPUTS . 

ENTER  INPUTS 
1*1,1, 1*1*1*1, 1,1*1 
DO  YO'J  1,1  ANT  TO  SEE  F  AND  F' 

1 =YES *  2=N0 . 

1 


•  0  0i2 1991 2269 1 648 


9  99498869 8 1 8 6 *  0 .  < .819141 2498 9 0 7 - 0 .  > 

0 9 9  * 498 8 6 9  8 1 8 6, 0 . >  < . 08191412498 9  0 7 ,  0 .  > 

9  98997  9897692 «  0.  :•  67  0992:2872725 ,  0. 

1  9979*59795 3 8 *  0.  >  < ,  i  841984774545.  0.  > 

9984972612258,  0.  :<  <'.54*6275487775*  0.  :• 

8  9  *54*  * i'i 8 8 677*  0 .  »  .  1  64 8 9  3  2 6 2 2: 2 2  2  <  0 .  :> 

9  97  9  9698  2562  8 *  0 .  ,  45 02 2 078 2777 6  •  0 .  1 

2991 *87924253  *  0.  :<  1  .  1  8009221  251  1  •  0.  - 

9  97 4*6:8566544 «  0.  >  1  .  26 28 026 069629 ,  0.  :> 

*  9974  84*8 2272,  0.  1  •. .  1  944  01  202481 4 *  0.  , 

9 96996980 2749 «  0 .  ■'<  2 0 2 1  014 2 8 4  2 0 5 *  0 . 

5  9  2 19  21  88 2  2  5  *  0 .  >  <: .  181 2  6  0  2: 5  7  0 5  8  2  *  0 . 

9  964972545*88,  0.  1  •  .  2 474627*  i  678 2*  0.  '< 

8  975481  4681  92*  0.  :«  <  .  1  722246193  748*  0.  ':< 

9 959*7*7*2 0 05 «  u.  >  •  .  202707758661 , 0.  :• 
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